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MIXING OF SUPERSONIC JETS INCLUDING THE EFFECTS OF TRANSVERSE 
PRESSURE GRADIENT USING DIFFERENCE METHODS 
byAnatole P. Kurkov 
Lewis Research Center 

SUMMARY 

The usual boundary- layer equations describing the steady- state mixing of parallel 
jets are supplemented by the momentum equation in the direction normal to the flow. 
This allows detailed computation of the flow field in the mixing region and simultaneous 
computation of the outer inviscid flow. 

An explicit and an implicit finite difference schemes have been developed and ap- 
plied in several illustrative examples. The examples include mixing of planar and axi- 
symmetric supersonic jets of different composition with both matched and unmatched 
static pressures. Numerical results were compared with available experimental data 
obtained for the unmatched- pressure case. 


INTRODUCTION 

Traditionally, jet mixing problems are solved using the boundary- layer equations. 
As a result, the transverse pressure field in the mixing region is assumed to be uni- 
form. In the case of a supersonic combustor this assumption may not always be justi- 
fied. Significant pressure gradients can be induced as a result of combustion, or they 
can be created at the injection port if the static pressures of the fuel and air are not 
matched. In such cases the transverse momentum equations must be considered. 

Previously (refs. 1 and 2), a numerical solution that incorporates the effects of the 
transverse pressure gradient was achieved by splitting the governing equations into a 
set of essentially inviscid equations and a set of viscous boundary- layer equations. In 
the first set of equations, which was written in the characteristic coordinates, viscous 
terms were treated as a perturbation. At each point it was necessary to iterate between 
the two sets of equations to match the two flows. 



In the present computational scheme, all equations are differenced in the same 
coordinate system and the variables associated with viscous and inviscid effects are de- 
termined simultaneously. Two finite- difference schemes, an explicit and an implicit, 
were developed and evaluated. 

The computer program was written for the solution of free or confined supersonic 
jet mixing. Initially, at the origin, the jets were assumed to have uniform properties 
and the turbulent Prandtl and Schmidt numbers were assumed to be unity. Figure 1 
shows the general configuration of the jets; geometry can be either planar or axisym- 
metric. 


Y 



Computations have been carried out for planar free jet mixing of hydrogen at Mach 
1. 67, and air at Mach 2.5. Both, matched and unmatched pressure cases were con- 
sidered. Additional cases were computed for the purpose of comparison with the theor- 
etical and experimental results from the literature. 

The computational scheme described in this report represents an alternate method 
to the usual scheme in which the inviscid and viscous flow fields are computed separ- 
ately. However, in case of reacting jets, the later method can not predict accurately 
the effects of combustion on the flow field. Schemes such as the one described in this 
report must then be used. In case of purely inviscid flows the present finite- difference 
scheme, in comparison with the method of characteristics, is easier to apply to the 
computation of ducted flows that involve multiple reflections from the boundaries. 
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SYMBOLS 


matrix whose element are a- t (eq. (31)) 

coefficient in eq. (24) 

matrix element defined by eq. (32) 

vector whose elements are bj (eq. (31)) 

coefficient in eq. (24) 

vector element defined by eq. (35) 

jet widths used in eddy viscosity expression, eqs. (36) and (37) 

coefficient in eq. (24) 

specific heat at constant pressure 

speed of sound 

defined by eq. (34) 

defined by eq. (33) 

defined by eqs. (29) 

total enthalpy 

static enthalpy 

J = 0 for planar, and J = 1 for axisymmetric geometry 
defined by eq. (21) 

Mach number 

coordinate normal to streamline 

grid spacing in n direction 

pressure 

velocity 

gas constant 

universal gas constant 

streamline coordinate 

grid spacing in s direction 

temperature 

velocity component normal to the direction of the shock 


Uj defined by eqs. (29) 

Vj defined by eq. (16) 

W k molecular weight for specie k 

X, Y rectangular Cartesian coordinates for planar configuration; axial and radial 

coordinates for axisymmetric configuration (see fig. 1) 

Yj half- jet height or jet radius at origin 

mass fraction for specie k 

/3 shock angle 

y specific heat ratio 

5 boundary- layer thickness 

e eddy (kinematic) viscosity; Mach angle, in the appendix 

Q flow angle 

M turbulent viscosity 

p density 

Subscripts: 

C centerline 

E external stream 

J central jet 

k refers to the particular specie; k = 1 for H 2 , k = 2 for Og, k = 3 for N 2 


ANALYSIS 

Basic Equations 

In the equations, the account of the turbulence is accomplished by using eddy vis- 
cosity, which is assumed to be a function of local flow variables and local jet param- 
eters. The variables in the equations are assumed to be time average quantities. 

Starting with full equations of motion in streamline coordinates and proceeding with 
boundary- layer simplifications, it is possible to derive 

pqjfl = -ig + J-/ M jsi) + j-^.. c -° se ia a 

3s 3s 3n \ 3n / Y 3n 
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as 9n 3 3n \ 3n/ 3n 3n 3s \ 3n/ 3 3n \ 3s/ 



where only the zero- and first-order terms were retained. In the equations (1) and (2), 
the turbulent viscosity is denoted by p; P, q, 0, and p denote pressure, velocity, 
flow angle measured from X axis, and density; s and n are streamline and the 
normal- to- the- stream line coordinates; Y is the radial distance from the axis; and 
J = 0 or 1, depending on whether the flow is planar or axisymmetric. The first equa- 
tion is the same as the usual boundary- layer equation. Considering the terms in this 
equation to be of zero order, all terms in the second equation are of the first order. 
Simultaneous solution of these equations would yield the solution accurate to the first 
order. However, the viscous terms in the second equation were found to affect the 
numerical solution in both the inviscid flow field and the mixing region only to a small 
degree. Omitting, therefore, the viscous terms in the second equation reduces it to 


2 30 3P n 

pq — + — = 0 

3s 3n 


(3) 


Although the simultaneous solution of the equations (1) and (3) throughout the flow 
field does not produce now a uniformly valid first-order solution, it correctly predicts 
the inviscid part of the flow. Consequently, the viscous region does not have to be nec- 
essarily small compared with the inviscid region; also the problem of matching of the 
two regions is eliminated. 

Assuming that the turbulent Prandtl and Schmidt numbers are equal to one, the 
energy and specie conservation equations are 

3H_ 1 3 / t 3H\ | J cos 0 3H ^ 

3s pq 3n V 3n/ pq Y 3n 


3s pq 3n \ 3n / pq Y 3n 


k = 1, ..., N 


(5) 


where H is the total enthalpy, H + q^/2, and and are the specie 

enthalpy and mass fraction. The continuity equation is given by 
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( 6 ) 


q — + p -^3 + pq — + J -^3 sin 6=0 
ds 3s 3n Y 


Denoting the temperature by T, universal gas constant by Rq, and the specie 
molecular weight by W k , the equation of state is given by 


P = pTR 0 


\ N a 

o) Z, 


^ w k 


(V 


It is possible to combine (ref. 3) the momentum equation (1), the continuity equation (6), 
and the equation of state (7) to obtain 


-UJC- i| — + 

pq 2 \c 2 I** 


3n CpT 1 9s / 7 K 3s \ q / pq L3n \ 3n/ Y 3nJJ 


.foyj. 
R Z^ w k 


da 


k _ j sin 9 


( 8 ) 


3s 


where the specific heat C , the gas constant R, and the velocity of sound c are defined 

ir 

by 




(9) 


■■'£3 


( 10 ) 


c V /2 

c =| 2_ rt 

IC^-R J 


(ID 
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Ignoring for the moment the terms on the right side of equation (8), it is seen that 
when the flow is supersonic, this equation together with equation (3) constitute a typical 
hyperbolic system. The terms on the right in equation (8), with the exception of the last, 
arise because of viscous diffusion. These terms originate from purely parabolic con- 
servation equations (1), (4), and (5). It is, therefore, seen that the equations are of 
mixed hyperbolic- parabolic character. 


The Explicit Method 

The difference equations for the parabolic equations (4) and (5) are obtained in the 
usual way. Written for the point (i, j) corresponding to the i^ interval in s direction 

i-L, 

and j Ln interval in the n direction they are 


Wj. 1 


^ . h u- h m-A 


As i,j piAj^ij-iM An ij 


An i,M 


K) - K) 

v H+u v k 4,j _ 1 


**i 


As i,j 


P i,j q i,j An i,3-l/2 


T u. cos 6 . . H. . - - H. . 
i J 11 1,3 1,3+1 i,]-l 

p. . q . - Y. . An. - + An. . .. 
ly J 1? 3 3 ^ ? j i > j“ 1 

Hm'Hi Hr K 1-1 


An i,j 


An. 


i, 3-1 


cos e t. i Mj±T Kh 


i Y i,i 




(12) 


(13) 


where As. . and An, . denote the step sizes in the streamline and the normal direc- 
1 \ J 

tions (fig. 2). In these equations viscosity is assumed to depend only on streamline 
coordinate. 

For the derivatives of P and 6 associated with the hyperbolic- type left sides of 
equations (8) and (3) a typical hyperbolic difference scheme (ref. 4, p. 262) is used: 



P i+l,j ~ P i,,j , e i,,j+l/2~ g i,j-l/2 = 


As- 


1,3 


An. 


3-1/2 


V 3 


(14) 
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where H’ and a ' denote differenced streamline derivatives of H and a given by 
equations (12) and (13). It is seen that 8 is computed from equations (14) and (15) 
midway between the grid points; other variables are computed at the grid points. For 
the purpose of the computation of coefficients in the difference equations, values of 8 
at the grid points, as well as the values of other variables between the grid points, are 
obtained by linear interpolation. 

In the axisymmetric case the two terms in equation (16) involving division by y 
become on the axis 


M i cose i.i < »i.wi-qi.)-i /Wi- %,i 


Y i,j in i,i + An i,J-l in i,i-l/2 \ a ”l,i 


An. 


(17) 


i, i- 1 


sin 8ii 8. • 


Ij p i. i-t-1/2 ~ e l,i-l/2 


Y l,i 


An. 


(18) 


i, 3- 1/2 


The last terms in equations (12) and (13) transform on the axis similarly, with q being 
replaced, respectively, by H and a. On the axis, also, normal derivatives must 
vanish so that 


H i+l,j-l - ^l.j+l 

K) = K) 

' K A+i, j-i ' k/ i+i,j+i 
q i+i,j-i =q i+i,j+i 


(19) 


and, since the flow angle is zero, 

0 i+l,j+l/2 0 i+l,j-l/2 (20) 

Conditions (19) and (20) also hold for the center plane in the planar case. 

The outer boundary of the jet is assumed to be reached when the deviation from the 
uniform external flow is sufficiently small. As the jet expands, additional points in the 
normal direction are included in the grid. The outermost streamline is assumed to be 
undisturbed flow. At each step, the equations are solved in the following sequence: (12), 
(13), (16), (14), and (15). After each advance in the streamline direction, the As and 
An are incremented using simple geometry. For example, 


9 


An i + l,j-l/2 ' An l,j-l/2 + As i, j(«i,J + l/2 - 6 i, j-l/ 2 ) 


The equations being explicit there are two necessary stability criteria to be satisfied. 
One is associated with parabolic equations and one with the inviscid hyperbolic parts of 
equations (14) and (15) (ref. 4): 


As i>j < K 


Pi i^i i 9 1 

1)] An. • K = — 


Mi 




1 

2 


As. ■ <J M 2 • - 1 An. • 
1,3 ? i,J 1 ,] 


( 21 ) 


( 22 ) 


where M is Mach number. In practice, in some cases, criterion (21) was not sufficient. 
A more stringent requirement was then imposed, 1/4 > K ^ 1/8. 


The Implicit Method 

The difference equations are transformed into implicit form by expressing the 
derivatives in the normal direction in terms of variables evaluated at the advanced dis- 
tance s. Thus, equation (12) becomes 

H i+l.i ' **t,j 1 **1 /%n.i+l ~ nj+i.i Hj+ij ' nj+i.i-i j 

As i,j P i, j q i, j An i, j-1/2 \ A “i,i-1 / 


J M l COa9 i,i ^1.1-1 




A) 


An i,i + An i,i-1 


(23) 


where, as in the explicit method, p is assumed to depend only on the streamline coor- 
dinate. This leads to a tridiagonal system of equations: 


- A, 


j H i+l> 3 


1 + B j H i+l,j " C j H i+l,j+l 


P i, 


j q i,j H i, j 


where 


A i= p i As i,i 


cos 6. . 
J 


An i,M/2 An i,j-l 


Y i,jKj-l + An i,3) 


(24) 


(24a) 
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°j = '-‘i as i,i 


1 , cosf, l..i 

An i, i- 1/2 a "i, j Y i,i( An i,i-l + an i,j) 


B. = p. -q. . + A. + C- 
3 f i,3 h i,3 3 3 

Another tridiagonal system results from the equations for 

‘ A i ( a k) + B i (“k) " C i f a k) = p i i q i i (®kl 

3 ' k 'i+l,j-l *'i+l,j n K/ i+l,j+l 1,3 1,3 ' K 'i, j 


(24b) 


(24c) 


(25) 


For the solution of systems of equations (24) and (25) there exists a simple algorithm 
(ref. 5, p. 14). 

Equation (14) now involves advanced values of 6 and q: 


1 f q i, j _ i\ P i+U ~ P i,.i . ^i+li 3+1/2 ~ p i+l, j- 1/2 = v 

2 \_2 ) As. , An 4 j 3 


- 6 . 


P i,3 q ij\ C l,3 


1,3 


3“ 1/2 


(26) 


V, = 



q i,3 + 


(S/jAi 

1 

q i,3 

P i,3 q i,3 


Mi 

J^i+lJ+l " VlJ Vi,) " q i+l,j-ll 


in i,i-l/2 

( an i,j a "i,j-l ) 

' Y M an i,j + an i,i-l j| 

> 



t sin e i»i 


W,. 


Y i,J 


(27) 


while equation (15) is unchanged. The differenced derivatives in the streamline direc- 
tion H' and in the identity (27) are calculated from the solution of systems of 
equations (24) and (25). 

It is possible to eliminate 6 from equations (15) and (26) to obtain 


E P. H ■ 

3 i+l,3- 


F S P i + l,j 


G j P i+l,j+l 


= U j - KL + l/2 V j in i, j-1/2 


(28) 
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where 


U 


rWu 


5+1/2 


(e. . , - e. )- Ml ~ j An i..i-i/2 p. 

\ i>j+l/2 ij 5 — 1/2/ / 2\ As i>l 

fa hi 


AS i,J 


- - KU 

1 Mi,|- 


i,j+l/2 As i,j-l/2 

An. . 1 
1/2 1}i ~ 1 


F. = 
] 


iS i,i + l/2 , 

K 

>i,j+l/2 As i, j-1/2 , ( M i,i ' 

£ 


ij 5+1/2 An i, j-1/2 

an i,J I 

K 

)- 1/2 1 

fa 

2 Y ■ 

h,j 

AS i,i _ 


G _ As i, 5+1/2 
3_ A “i,i 


/ (29) 


The unknowns in equations (28) are P and q (which appears in Vj); however, P can 
be eliminated using the implicitly differenced form of momentum equation (1), 

P; . , , - P, ; 


n a - - IhhlLlhl 

j As As An 

AS i,j AS i,j An i> j- 1/2 


+ 


i q i+l j j+1 ~ q i+l,.i _ q i+l,j ~ q i+l,j-l ) 


an i,i 


in ij-i 


J “l cos ^i,j q i+l, j+l ~ ^i+l. j-l 


Y U 


in i,i-l + in i,j 


(30) 


The result is a five- diagonal system of equations. In matrix notation it is 

Aq = B (31) 


The elements of matrix A, a- v are zero except for j - 2 < k ^ j + 2: 

J,K 
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a i,j-2 - E J A j-l 




a. . = E-C. h - D-B. + G.A. . - d. L 

3,3 3 3- 1 3 3 3 3+1 3 f 


a* • 4 = D‘C- - B. i G* 
3, 3+1 3 3 3+1 3 


a j , 3+2 - C j+l G j 


(32) 


where is defined by 


, A. 1+1/2 An i,i-l/2 q ‘’i + fcp\.i Tl ’i 


3 


AS i,3 




(33) 


and 


d i 

D. = F. 1 

3 3 


p i,3 q i,3 


(34) 


The general term for vector B, bj is given by 


bj = u r (pq 2 ) i>j+ i /2 An i, 3-1/2^77—^ 


HiM 

1,3 

R oV' 

(“A, 


w k 

k 


f j[ p m + 

(pq 2 ) 




(35) 
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Denoting the lower boundary streamline by jp and the outer boundary streamline by j^, 

it can be seen from equations (32) and (35) that elements a- . and b- must be modi- 

J y -K J 

fied for j ^ jp + 2 and j > - 2 in order to satisfy proper boundary conditions. At 

the outer boundary they are evaluated assuming that the flow variables for j > are 
equal to the free- stream values. The is chosen so that the flow disturbance relative 
to the free stream is small. 

The treatment of the lower boundary is the same for distances X that are less than 
the potential core length. Beyond the potential core, identities (17) and (18) must be 
used on the axis, with q evaluated at the advanced level in (17). In both, axisymmetric 
and planar cases, the matrix coefficients near the center streamline are evaluated using 
equations (19) and (20). The system of equations (31) was solved using specialized Gauss 
elimination method. 


RESULTS AND DISCUSSION 
Illustrative Examples 

The results were obtained for free, two-dimensional jet mixing with both matched 

and unmatched pressures and for two different turbulent viscosity models. At the origin 

both jets were assumed to have uniform properties. The center jet was hydrogen with 

P = 1. 01x10 or 3. 03X10 newtons per square meter, M = 1. 67, q = 2220 meters per 

5 

second, and T = 306 K. The external stream was air with P = 1. 01x10 newtons per 
square meter, M = 2. 48, q = 1630 meters per second, and T = 1110 K. The two turbu- 
lent viscosities considered were constant viscosity and a more realistic formulation 
based on Prandtl kinematic eddy viscosity model (ref. 6, pp. 599 and 607): 

e = 0. 014 bp j q E - qj + 0. 0014 (36) 

for X less than the length of potential core, and 

e = 0. 037 bjy 2 qg - q c (37) 


for X greater than the potential core length. In the equation (36) bp j is the width of 


= 0. 1 and 


the mixing zone measured between points where | a. j - = 0. 1 and | - a gp| = 

0. 1. The mass fraction a ^ designates the component originally constituting the jet 
(hydrogen). This definition of b^ j differs somewhat from one given in reference 6 
where bp ^ is defined using the velocity difference. Mass- fraction difference was cho- 
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sen here because the change in composition is clearly related to viscous effects, whereas 
the change in velocity may be due to inviscid effects. The width b^^ is measured be- 
tween the jet centerline and the point where q - q^-. | =1/2 j - q^, ( . 

The value chosen for the constant viscosity case was p = 4. 79X10” ^ newtons- 
seconds per square meter, which is of the order of magnitude found in the fully devel- 
oped turbulent jets. 

The results were obtained for three cases. In case I the pressures were matched 
and the viscosity was constant. Case II was calculated with hydrogen pressure three 
times the external stream pressure while, as in the first case, viscosity was constant. 

In case HI pressures were matched and eddy viscosity was calculated from equations (36) 
and (37). 

Figure 3 presents pressure profiles at downstream distances X/Yj = 0.2 (fig. 3(a)); 

Implicit, fine grid 

Implicit, coarse grid 

a Explicit, fine grid 
o Explicit, coarse grid 



(b) Downstream distance, X/Yj, 1.0 and 3.0. (c) Downstream distance, X/Yj, 40 and 100. 

Figure 3. - Pressure profile for constant viscosity. Initial pressure ratio, PyP E , 1. 
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X/Yj - 1 and 3 (fig. 3(b)); X/Yj = 40 and 100 (fig. 3(c)). The half- jet height Yj was 
assumed to be 0. 1905 centimeter. Since the pressures were matched initially, the pres- 
sure waves noted in these figures are due to mixing. Close to the injector (fig. 3(a)), 
the pressure waves are steep; however, further downstream the wave amplitudes are 
progressively attenuated (figs. 3(b) and (c)). 

In figure 3 comparison was also made between implicit and explicit solutions, with 
the degree of convergence indicated in figures 3(a) and (b) by presenting results com- 
puted for two grid sizes. It is seen that the convergence of the implicit method is good 
except in the narrow regions around the pressure peaks. The fine grid spacing was half 
the value of the coarse grid spacing in both implicit and explicit solutions. Decreasing 
the grid spacing further was considered impractical. The actual mesh size was varied 
during the calculation; however, past the initial portion, the number of points in the 
normal direction was limited to 90 for coarse grid and to 180 for fine grid. The explicit 
solution at X/Yj = 0. 2 (fig. 3(a)) is quite irregular with oscillatory behavior in the 
hydrogen region (Y/Yj < 1). In this region no particular improvement could be noted by 
using a finer grid. However, further downstream, as the disturbed region propagates, 
the explicit solution approaches a smooth curve. Using this method, the resolution of 
the steep pressure wave in the air region is furthermore seen to be superior to the reso- 
lution obtained using the implicit method. At a distance of X/Yj = 100 (fig. 3(c)) the 
difference between the explicit and implicit solutions appears to be significant on a rela- 
tive basis although the absolute magnitude of the pressure wave at this distance is small. 
The agreement in the wave speed, however, is very good. 

Another indication of the accuracy of the implicit method is presented in the next 
section where it is applied to the calculation of the inviscid flow field and compared with 
the method of characteristics. 

The oscillatory behavior in the explicit solution was observed whenever criteria (21) 
required that As « An. This is usually the case in the region close to the injector exit 
where An must be small so that steep pressure waves could be resolved. The oscilla- 
tion was also observed in the implicit solution when As « An. In this case, however, 
there is no restriction on As so that in figure 3 and in the following figures As was 
chosen equal to An. 

Figure 4 presents profiles of hydrogen mass fraction for several downstream dis- 
tances. It is seen that the implicit solution follows the explicit closely and that there is 
only a small difference between these and the solution obtained using conventional 
boundary- layer equations (P held constant). It is noted that the mixing zone is consid- 
erably narrower than the width of the pressure field. 

The results for case n, when the jet pressure is three times that of the external 
flow, are shown in figures 5 and 6. The mixing region assumes a pressure somewhat 
below the average, with the expansion wave moving in the direction of the central jet and 
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Dimensionless vertical distance, Y/Yj 


Figure 4. - Profiles of hydrogen mass fraction constant 
viscosity. Initial pressure ratio, 1. 


Downstream 



0 .5 1.0 1.5 2.0 2.5 3.0 


Dimensionless vertical distance, Y/Yj 

Figure 5. - Pressure profiles. Constant viscosity. Initial 
pressure ratio, Pj/P E , 3. 
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Figure 6. - Profiles of hydrogen mass fraction. Constant viscosity. 
Initial pressure ratio, Pj/P^ 3. 


the compression wave moving in the opposite direction. In figure 6 the mixing region is 
displaced considerably in the direction of the lower pressure external flow. 

Figures 7 and 8 were obtained using equations (36) and (37) for the eddy viscosity. 
Physically, this represents a more realistic assumption and is typical of a number of 
available models that attempt to correlate experimental results. The amplitude of the 
pressure waves is much smaller in this case, although, a mild and fairly uniform pres- 
sure increase can be noted at X/Yj = 3 . It is evident that the near pressure field is 
very strongly dependent on the magnitude of the viscosity. 


Downstream 

distance, 

X/Yj 



o 

rj 


o I I I A i 1 • 

' 0 .5 1.0 1.5 2.0 2.5 3.0 

Dimensionless vertical distance, Y/Yj 


Figure 7. - Pressure profiles. Prandtl viscosity. Initial 
pressure ratio, Pj/P E , 1. 
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0 .5 1.0 1.5 2.0 2.5 3.0 

Dimensionless vertical distance, Y/Yj 


Figure 8. - Profiles of hydrogen mass fraction. Prandtl 
viscosity. Initial pressure ratio, Pj/P^, 1. 



Figure 9. - Hydrogen mass-fraction distribution along center 
plane. 
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Figure 9 shows the distribution of hydrogen mass fraction along the center plane for 
all three cases. The potential core length for the Prandtl viscosity case is greater than 
that for the constant viscosity case due to smaller viscosity close to the injection point. 
The difference in the potential core length and in the rate of hydrogen concentration decay 
for the unbalanced pressure case is mainly due to the expansion of the hydrogen jet. 

Figure 10 illustrates the positions of the mixing zones for cases I and n as deter- 
mined by 1 and 99 percent hydrogen boundaries. Relative to the matched pressure case, 
the mixing zone in the unmatched pressure case is displaced in the direction of the lower 
pressure external flow. 

initial 

pressure 



Figure 10. - 1 And 99 percent jet boundaries. 

Comparison With the Literature 

In figures 11 to 14 results were obtained for several axisymmetric jet mixing prob- 
lems for the purpose of comparison with the numerical method of reference 1. Figure 11 
presents pressure distribution along the centerline for free jet mixing assuming matched 
pressures and a constant turbulent viscosity. The present solution (explicit method) 
produces a much steeper pressure peak on the axis. A similar conclusion is reached 
from figure 12, which was obtained for underexpanded jet exhausting in a ducted coaxial 
flow assuming purely inviscid interaction. In this figure, the pressure peak at about 
X/Yj = 18 is due to the reflection of the original compression wave from the upper wall. 
Included in the figure is the implicit solution computed for two grid sizes and also the 
solution obtained using the method of inviscid rotational characteristics presented in the 
appendix. It appears that a very fine mesh spacing is required for an accurate resolution 
of the pressure waves on the axis. There is still a noticeable departure of the fine- grid 
implicit solution from the points plotted for characteristics solution around the minimum 
pressure point in the first expansion wave. However, the difference between the char- 
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Figure 11. - Coaxial air jets, axial pressure distribution. Viscosity, 
2.49xl0" 2 newton-second per square meter,- jet radius, 2.5 centi- 
meters. Initial properties of jet: pressure, 1. 013xl0 5 newtons per 
square meter; Mach 2; temperature, 1000 K. Initial properties of 
external flow: pressure, 1.013xl0 5 newtons per square meter; 
Mach 3; temperature, 300 K. 



Figure 12. - Underexpanded ducted jet, centerline pressure distribu- 
tion, inviscid case. Jet conditions: Mach 2.0; pressure, 5. 06xl0 5 
newtons per square meter; temperature, 1100 K; gas, hygrogen; 
jet radius, 0.1 centimeter. External flow conditions: Mach 3.38; 
pressure, 2.533 newtons per square meter; temperature, 1500 K; 
gas, air; duct radius, 0.376 centimeter. 
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Figure 13. - Underexpanded ducted jet, pressure profile. 
Downstream distance, X/Yj, 4.0; inviscid case. 



Figure 14. - Initial region of underexpanded ducted jet. 
Centerline pressure distribution; viscosity, 2. 15x10 
newton-second per square meter. 


acteristics solution and the results obtained in reference 1 is even more pronounced. 
Unfortunately, no information is given in this reference on the degree of convergence of 
the solution. 

A further indication of the accuracy of the implicit solution can be obtained from 
figure 13 where the pressure distribution at X/Yj = 4. 0 is compared with the corre- 
sponding results computed using the method of characteristics. It is seen that implicit 
method tends to smooth out the discontinuous compression wave at the outer boundary of 
the disturbed region. Similar tendency was noticed earlier in the discussion where im- 
plicit and explicit methods were compared. 

Additional comparison with results from reference 1 is made in figure 14 which pre- 
sents the centerline pressure distribution in the initial region of underexpanded jet in a 
ducted flow assuming a constant viscosity. It is difficult to explain the pressure rise at 
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X/Yj = 1. 6, noted in the solution of reference 1 even when considering (in view of results 
in fig. 12) that the grid spacing for this case was perhaps too coarse. Despite such poor 
agreement with reference 1 in figure 14, relatively small difference (about 10 percent 
for the range of X/Yj covered in the figure) was observed comparing the hydrogen 
mass- fraction distributions on the centerline. 


Comparison With the Experiment 

In this section, the experimental results from a recent study of wall slot injection 
into a supersonic air stream (ref. 7) are compared with the finite- difference solution. 

A sketch of the wall slot is presented in figure 15. Also shown in this figure is the gen- 


Free 

stream 



Figure 15. - Wall slot jet. Slot height, 1.27 centimeters; lip thickness, 0.0127 centimeters. 


eral wave pattern and the extent of the boundary layer above the splitter plate for the 
particular series of experiments. The bottom-wall boundary layer was small compared 
with the slot height, and it did not seem to affect appreciably the wave pattern. In these 
experiments the composition of the jet was air, as was the composition of the free 
stream. The jet was supersonic with M = 1. 98 and total temperature 297 K. In the 
free stream the total pressure and temperature were ambient with M = 4. 19. 

In the finite difference solution, viscosity was assumed to be given by 

P = pe + P’ 

5 

where p' was taken to be of the order of molecular viscosity, p* = 0.744X10 newtons- 
seconds per square meter. The eddy viscosity e was assumed to be given by equations 
(36) and (37) with the constant in equation (36) deleted and bg ^ increased by adding a 
constant equal to the boundary- layer thickness. The computation was also carried out 
taking e = 0 and, consequently, p = p’; however, the difference was very small. The 
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(a) Initial pressure ratio, Pj/P E , 1.31. 



(b) Initial pressure ratio, Pj/P E , 2.5. 

Figure 16. - Pressure distribution on bottom wall, initially uniform 
external flow. 


boundary layer on the bottom wall was not included in the calculation, and in figures 16 
and 17 neither was the boundary layer above the splitter plate. 

Figure 16(a) presents the wall- pressure distribution for the under- expanded slot jet 
with Pj/Pg = 1.31 at the injection point. For X/Yj < 4. 0, the agreement with the ex- 
periment is not very good; however, further downstream the agreement is much better. 
At X/Yj about 5. 5, in both the theoretical and experimental results, there is evidence 
of a weak reflection of the expansion wave from the mixing zone . 

Figure 16(b) presents the wall- static pressure distribution for a higher initial pres- 
sure ratio of the two jets, Pj/Pg =2.5. The contours of the dividing streamline for this 
case are shown in figure 17 where the experimental points were estimated from the 
photographs of reference 7. The agreement with the experiment is somewhat better in 
this case. 

In order to explain the over- expansion and the higher apparent wave speed observed 
in the experimental results in figure 16(a), the computation was carried out with the eddy 
viscosity increased considerably beyond the value originally assumed. With this large 
viscosity it was possible to match the experimental wave speed and to obtain qualitatively 
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Figure 17. - Contour of dividing streamline, initially uniform exter- 
nal flow. Initial pressure ratio, Pj/P E , 2.5. 



Figure 18. - Assumed velocity distribution in bound- 
ary layer. 


similar pressure distribution. However, the mixing region in this case was much wider 
than indicated in the schlieren photographs of reference 7. 

Next, the effect of the boundary layer above the splitter plate was investigated. It 
was assumed that density varies linearly in the boundary layer as given by 


P =P E °* 2 + 



which approximately matches the density variation in figure 44 of reference 7. To avoid 
the transonic and the subsonic regions in the boundary layer, uniform flow was assumed 
between the wall and the point where M =1.5. The total enthalpy was assumed to be 
uniform throughout the boundary layer. The resulting velocity distribution in the bound- 
ary layer is shown in figure 18. 

It can be seen in figure 19 that in this case the static pressure distribution more 
nearly approximates the experimental data. The apparent discrepancy in the wave speed 
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Figure 19. - Pressure distribution on bottom wall considering bound- 
ary layer in external flow. Initial pressure ratio, 1.31. 


that is still present is probably due to the upstream propagation of the pressure disturb- 
ance on the bottom wall through the subsonic portion of the boundary layer. 


CONCLUDING REMARKS 

A numerical procedure has been developed for simultaneous computation of viscous 
and inviscid flow fields resulting from the supersonic mixing of jets. 

The procedure can use either an explicit or an implicit finite difference scheme, 
each one having certain advantages and disadvantages. The explicit scheme sometimes 
produced an oscillatory solution. The oscillations were noticed in cases when the para- 
bolic stability criterion requires that the streamwise grid spacing be much smaller than 
the spacing in the normal direction. The implicit scheme is generally free of this be- 
havior; however, it is found that the resolution of the pressure waves generated as a 
result of mixing is not as good. The tendency is to smooth out the pressure gradients 
in the region where they are steep and, therefore, to reduce the magnitude of the pres- 
sure peaks. The mixing process, however, is found to be insensitive to the degree of 
accuracy with which the steep pressure variations are described. Illustrative examples 
comprise jet mixing with matched and unmatched static pressures. As a limiting case, 
the example for unmatched static pressures includes flow computation for purely inviscid 
jets. For this case, also, the method of inviscid rotational characteristics was used. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, September 29, 1971, 

764-75. 
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APPENDIX - METHOD OF CHARACTERISTICS FOR UNMATCHED-PRESSURE JETS 


The inviscid rotational characteristics equations are derived in terms of variables 
P and 9 from equations (3) and (8) in which the viscous terms are deleted. For rota- 
tional flows, working with these variables is more convenient than with usual variables 
q and 9 or Prandtl- Meyer function and 9, which necessitate the introduction of another 
variable such as entropy. 

Assuming that the data at mesh points 1 and 2 are known, the solution for point 3 is 
obtained from 


Y 3 - Y 1 = (X 3 - Xj) tan - e) 

Y 3 — *^"2 — (X3 — X2 ) (#2 ^2^ 

sin e, cos e. sin e, sin 0i(Xq - X-.) 

9 , - 9 q + (Pq - P, ) + J — = 0 

1 J yjPj cos (0J - e J )Y 1 


sin e 2 cos £2 
^ 2 P 2 


(Pq ~ P 9 ) + J 


sin €2 s i n ® 2^3 “ 
cos (0g ^ 2 )Y 2 


= 0 


(Al) 

(A2) 

(A3) 


(A4) 


where e denotes the Mach angle and y the specific heat ratio. 

In the example in figures 12 and 13 the expansion of the central jet and the compres- 
sion of the outer flow are determined so that the respective final pressures and final 
flow angles are the same. Since the gas is assumed to be calorically imperfect, the 
Prandtl- Meyer function for the expansion wave cannot be expressed as a function of 
Mach number only. It is more direct to solve the expansion wave by integrating 9 first, 
since it can be expressed as a function of T only: 


d 9 



2[H - h(T)] 
y(T)RT 


1 


c p ( t) 
2[H- h(T)] 


dT 


(A5) 


where for h(T) and Cp(T) the same expressions were used as in the finite- difference 
program. In the expansion region pressure can also be expressed as a function of T 
by integrating 

jT, CJT) 

^=— E dT (A6) 

P RT 
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For the compression shock in the external stream it can be shown, using the con- 
servation of mass and momentum in the direction normal to the shock and the equation 
of state, that the following equation must hold: 


T 


1 




+ ^|K 

+ r /u 2 


(A7) 


where 1 and 2 refer to the upstream and downstream sides of the shock and U de- 
notes the component of q normal to the direction of the shock. Denoting the shock 
angle by /3, and U 2 are given by 


U 


1 " 



sin /3 


(A8) 


and 


U 


2 _ 



(A9) 


In the computation scheme, initially /3 must be assumed. Using this value of fi, 

Tg is then iterated until equation (A7) is satisfied within acceptable limits. This allows 
calculation of 0 2 from 


sin (/3 - 



the density p 2 from the continuity equation across the shock, 

Pl u l = p 2 u 2 


(A10) 


(All) 


and P 2 from the equation of state, 


P 2 = p 2 RT 2 (A12) 

In the expansion region, the integration of equation (A5) gives 0 = 0(T) in tabular 
form. The temperature at the end of expansion is obtained by interpolation from this 
table using 0 = 0 2 . Integration of equation (A6) gives now the pressure at the end of 
expansion. This pressure is compared with P 2 . If the difference exceeds a prescribed 
tolerance, a new value of /3 is assumed, and the iteration is continued. 
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Downstream of the initial point the treatment of the shock differs to the extent that 
/3 is iterated until the compatibility condition (A4) is satisfied using 6 and P deter- 
mined from equations (A10) and (A12). 

At the interface of the two jets, the pressure and flow angle must be the same on 
both sides, and they must also satisfy the compatibility equations: on the side of the 
outer stream, equation (A3); and on the jet side, equation (A4). 
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